Band convergence and linearization-error correction of all-electron GW calculations: 

the extreme case of zinc oxide 
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Recently, Shih et al. [Phys. Rev. Lett. 105, 146401 (2010)] published a theoretical band gap for 
wurtzite ZnO, calculated with the non-selfconsistent GW approximation, that agreed surprisingly 
well with experiment while deviating strongly from previous studies. They showed that a very large 
number of empty bands is necessary to converge the gap. We reexamine the GW calculation with 
the full-potential linearized augmented-plane-wave method and find that even with 3000 bands the 
band gap is not completely converged. A hyperbolical fit is used to extrapolate to infinite bands. 
Furthermore, we eliminate the linearization error for high- lying states with local orbitals. In fact, 
our calculated band gap is considerably larger than in previous studies, but somewhat smaller than 
that of Shih et al,. 
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In a recent Letter Shih et al} presented a new cal- 
culation for the band gap of ZnO in the wurtzite struc- 
ture employing the GW approximation^ for the electronic 
self-energy. They used a conventional non-selfconsistent, 
one-shot approach^il in which neither the quasiparti- 
cle energies nor the quasiparticle wave functions are up- 
dated. Instead, the GW self-energy E'^(r,r';a;) [Eq. ^] 
is constructed from the Kohn-Sham Green function taken 
from a density-functional theory calculation with the 
local-density approximation (LDA) for the exchange- 
correlation energy functional. The quasiparticle energy 
^nq with band index n, Bloch vector q, and spin a is 
then obtained from the nonlinear equation 



Kq = <.q + (^^q|S^(Kq) ^ Cl^^^q) 



(1) 



with the Kohn-Sham energy ej^^, the wave function 
</7j^q(r), and the exchange-correlation potential v^^{r). 
Offdiagonal elements of E'^ — v^^ are neglected. 

All previous calculations^"— invoking the one-shot 
LDA+GW approach showed that the band gap of 
wurtzite ZnO is underestimated with respect to the ex- 
perimental value by more than 1 eV. They fall in the 
range 2.12-2.6 eV while the experimental gap amounts 
to 3.6 eV,^ after correction for vibrational effects. This 
large underestimation is untypical for GW calculations 
of sp-bound systems. 

The Letter of Shih et al. addressed two issues: first, 
the erroneous hybridization effects between Zn 3d and 
O 2p states that results from the self-interaction error 
within the LDA^ and second, the band convergence in 
the correlation part of the self-energy. The first problem 
was tackled with the LDA+J7 approachjii in which an 
orbital-dependent potential corrects the position of the 
3d bands and, thus, reduces hybridization effects with 
the O 2p states. However, the combinatioir LDA+J7 
and GW yields a band gap that is still well below the 
experimental value. Therefore, the authors investigated 
the secoird issue by carefully converging the correlation 
self-energy and the dielectric matrix with respect to the 
number of bands. They performed calculations with up 



to 3000 bands corresponding to a maximal band energy 
of 67 Ry as well as a cutoff for the dielectric matrix of 
up to 80 Ry and showed that the resulting GW band 
gaps, 3.4 eV for LDA+GVF and 3.6 eV for LDA+U+GW, 
turned out to be in very good agreement with experi- 
ment. They also demonstrated that a too small energy 
cutoff for the dielectric matrix can lead to a false con- 
vergence behavior: the band gap seems to converge with 
many fewer bands, but toward a value that is too small. 

These new results for ZnO are in striking contrast to 
previous studies. If they are correct, they cast doubt on 
all GW calculations published so far, not only for ZnO 
but also for other materials, especially for systems with 
localized states. In fact, Shih et al. point out at the end 
of their paper that "many of the previous quasiparticle 
calculations ... involving localized states may need to 
be reexamined." This will likely provoke a controversial 
debate that reaches different fields of solid state physics 
and requests a rapid clarification. 

Since its publication the paper has aroused criticism, 
mostly from the all-electron community who pointed out 
that calculations with the pseudopotential approxima- 
tion cannot give a definite answer for the real GW gap be- 
cause of the neglect of the core- valence exchange effects, 
the pseudized form of the valeirce wave functions, and 
the inaccurate description of high-lying states. There- 
fore, they attributed the large difference between the new 
result and the previous all-electron calculations^"^ (2.12- 
2.44 eV) to the approximations inherent to the pseudopo- 
tential approach or to numerical errors in the calculation. 
It is, thus, vitally important that the GW band gap of 
ZnO is reinvestigated and thoroughly coirverged with a 
genuine all-electroir method to provide a standard for the 
discussion that will follow. 

In this paper, we present an all-electron LDA+GM^ 
calculation for wurtzite ZnO that is based on the full- 
potential linearized augmeirted-plane-wave (FLAPW) 
method. For simplicity, we restrict ourselves to the 
standard LDA approach for the noninteracting starting 
point without an additional Hubbard U parameter. Our 
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calculation yields a band gap that is, in fact, much larger 
than that of the previous calculations, but somewhat 
smaller than the result of Shih et al. . We go beyond their 
approach in two respects: we employ neither the pseu- 
dopotential approximation nor a plasmon-pole model for 
the dielectric matrix. Instead, the screened interaction 



M^(r,r';w) = w(r,r')+ jj v{r,r")P{r" ,r"';uj) 

xW{r"\r';oj)d^r"d^r"' (2) 

is calculated explicitly within the random-phase approx- 
imation for the polarization function 

X G" ir\r, uj' -uj)duj' , (3) 

where the Green function G"^ (r, r', lj) is constructed from 
Kohn-Sham energies and wave functions. The frequency 
convolution of the self-energy 



E-(r,r';c.) = 



I 

2^ 



G''{r,r';u;+u;')W{r, r'; 



duj' 
(4) 

{•q is a positive infinitesimal) is evaluated analytically for 
w(r, r') [see Eq. and with a contour integration^ii^ 
on the complex frequency plane for the remainder 
W{y, r'; w) — w(r, r'). The nonlinear equation ([IJ is solved 
on an energy mesh around e^q with a cubic spline inter- 
polation between the mesh points. Thus, no additional 
Taylor expansion of the self-energy is needed. Details of 
the implementation can be found in Ref. IT^ . 

We carefully converged the number of empty bands 
for the calculation of both the polarization function and 
the correlation self-energy as well as the parameters 
for the all-electron mixed product basis,— in which 
we represent the dielectric matrix. While the ground- 
state electron density was converged with a standard 
LAPW basis with moderate cutoff parameters (^max — 8, 
Gniax ~ 4.3 flg"^, where ag is the Bohr radius), we had 
to employ much larger cutoffs to generate enough wave 
functions for the GW calculation: a linear momentum 
cutoff of Ginax = 8.0 Cj^^ and an angular momentum cut- 
off in the muffin-tin (MT) spheres of ^max = 12. Further- 
more, in order to avoid linearization errors in the MT 
part of the LAPW basis^i^iil we added local orbital a^^'^^ 
(LOs) with different angular momentum quantum num- 
bers and energy parameters distributed over the relevant 
energy range: 292 LOs for Zn (five LOs for each Im chan- 
nel with I = 0-3; three for I — 4. two for I = 5, and one 
for I — 6) and 186 for O (four LOs for / = 0-3, two for 
I = 4, and one for 1 — 5). We also treat the 3s and 3p 
semicore states of Zn explicitly with LOs. 

For the mixed product basis we found an angular mo- 
mentum cutoff of 4 in the spheres and a suprisingly small 
linear momentum cutoff of 2.4 Aq ^ in the interstitial re- 
gion to be suflacient. However, we had to take into ac- 
count many product functions in the MT spheres, which 
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Figure 1: Band convergence of the quasiparticle band gap 
of ZnO employing a 4x4x4 k-point set and calculated with 
(pluses) and without local orbitals (LOs) (crosses) for high- 
lying states. The solid lines show the hyperbolical fits. We 
also indicate results with finer k-point samplings (stars) cal- 
culated with LOs and 500 bands. The dashed lines show the 
hyperbolical fit shifted to align with these results. The fit 
asymptote for the 8x8x8 k-point set at 2.99 eV (dotted line) 
is considered the best estimate for the all-electron one-shot 
GW band gap. 



after optimization led to 177 MT functions for Zn (ten, 
eight, eight, seven, and six radial functions per Im chan- 
nel for I — 0-4, respectively) and 127 for O (eight, six, six, 
five, and four radial functions for I — 0-4, respectively). 
Obviously, the rapid variations close to the atomic nuclei 
must be accurately described. Within the all-electron 
mixed product basis this is possible with a relatively 
modest number of MT functions, while in a pure plane- 
wave approach a very large number of plane waves is 
necessary to resolve the variations adequately. This ex- 
plains the finding of Shih et al. that the dielectric matrix 
must be converged to very large energy cutoffs. The total 
number of mixed product basis functions in the calcula- 
tions is less than 700 per k point. This number is further 
reduced to around 490 by constructing linear combina- 
tions that are continuous in value and radial derivative 
at the MT sphere boundaries.^*' 

Figure [T] shows the quasiparticle band gap of ZnO as a 
function of the number of states included in the calcula- 
tion of the polarization function and the correlation self- 
energy. The calculations were performed with a 4x4x4 
k-point sampling of the Brillouin zone. There is a large 
difference between calculations with (pluses) and with- 
out the LOs for unoccupied states (crosses), which shows 
the importance of eliminating the linearization error of 
the conventional LAPW basis. As the linearization er- 
ror becomes larger for higher and higher bands, it is not 
surprising that the difference between the convergence 
curves grows toward increasing numbers of bands. We 
find an asymptotic difference of 0.5 eV. The calculations 
without LOs for unoccupied states already converge with 
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Figure 2: Representation error A(e) (see text) of the set of 
spherical MT functions at an oxygen atom with (solid line) 
and without local orbitals (LOs) (dashed line) as a function of 
energy. The Fermi energy is set to e = 0. The "spikes" clearly 
identify the positions of the energy parameters. In the case 
without LOs, the representation error soon approaches unity 
(dotted line) for energies above 10 Ry, which means that the 
basis becomes practically orthogonal to the exact solution of 
the radial scalar-relativistic Dirac equation. 



a few hundred bands, which could have led the authors 
of the previous all-electron studies to believe that their 
calculations are sufficiently converged. In fact, when the 
converged value of 2.27 eV is corrected with respect to 
finer k-point samplings, we arrive at 2.44 eV, which lies 
at the upper edge of the range of the all-electron GW 
band gaps published so far. 

The linear momentum cutoff allows the interstitial part 
of the LAPW basis to be converged in a systematic 
way. In the MT spheres, however, the basis is linearized 
around predefined energy parameters, which gives rise to 
the linearization error for high- lying states. We now an- 
alyze this error in more detail and show that it can be 
eliminated very effectively with the LOs. Figure [2] shows 
the representation error A(e) of the MT basis at an oxy- 
gen atom for the angular momentum I — as a function 
of energy. We define 



A(6) 



(5) 



where S is the oxygen MT sphere radius, R{e, r) is the 
normalized solution of the radial scalar-relativistic Dirac 
equation (cf. Ref. [13) for the energy e and the angular 
momentum I = 0, and i?r(e,'') is its best representation 
in terms of linear combination of the MT functions con- 
tained in the LAPW basis. In the conventional basis 
there are only two radial functions available: the solu- 
tion -R(epar,'') of the radial Dirac equation for / = at 
the energy parameter Cpar, in this case Cpar — —1.24 Ry, 
and its energy derivative (ii?(epar, f) / f^Epar- [In the litera- 
ture these functions are commonly denoted by u;=o('') 
and ui^Q{r), respectively.] The dashed line in Fig. [5] 



shows that this conventional basis represents the occu- 
pied states very accurately, which are located in the small 
energy range —1.33 < e < ORy, but fails to describe un- 
occupied states at higher energies. In fact, at energies 
above 10 Ry the exact solution i?(e, r) becomes practi- 
cally orthogonal to the MT basis, and the representa- 
tion error approaches unity. Interestingly, at this energy, 
which roughly corresponds to the 500th band, the con- 
vergence curve for the calculations without the LOs levels 
off, and the band gap seemingly converges (see Fig. [Ij . 
The description of wave functions at higher energies is 
considerably improved, if LOs are added to the LAPW 
basis. The solid line shows the corresponding represen- 
tation error with four additional LOs at 12, 24, 36, and 
60 Ry. As can be seen, the error remains below 10"'^ for 
a very large energy range up to 65 Ry. 

As Fig. [1] shows, the calculations with eliminated lin- 
earization error (pluses) take far more bands to converge. 
In fact, even with 3000 bands the band gap is still not 
completely converged. Therefore, we fitted the values 
with the hyperbolical function 



f{N) 



N ~Nf 







(6) 



where a, b, and iVo are fit parameters. It is surprising 
how closely the fitted curve (solid line) follows the calcu- 
lated data points. This makes us confident in taking the 
fit asymptote b as the band gap extrapolated to infinite 
bands. Furthermore, we have recalculated the band gap 
with finer k-point meshes (up to 8x8x8) and 500 bands 
(crosses). The dashed lines show correspondingly shifted 
hyperbolical fits. The asymptote of the fit corresponding 
to an 8x8x8 k-point sampling is found at 2.99 eV, which 
we take as the final best estimate for the all-electron one- 
shot GW band gap. 

This band gap is 0.4-0.9 eV larger than the previously 
reported values. Both the large number of bands that 
are needed for proper convergence and the elimination of 
the linearization error, which has not been undertaken in 
the previous all-electron studies, are responsible for this 
large difference. Our value is still about 0.4 eV smaller 
than the band gap of Shih et ai, though. In fact, a cer- 
tain discrepancy should be expected because of the pseu- 
dopotential approximation and the plasmon-pole model 
for the dielectric function used in Ref. [l]. The pseudopo- 
tential approximation not only neglects the important 
contribution of core- valence exchange. It also yields ac- 
curate wave functions only in the vicinity of the atomic 
electron energies, but not for high- lying states. This er- 
ror is very similar in spirit to the linearization error of the 
LAPW basis and is also of the same mag nitude,ii Thus, 
it should have an impact on the GW results comparable 
in size to the linearization error, whose elimination gives 
rise to a sizable correction of 0.5 eV, as we have shown 
in the present work. 

With the LDA band gap of only 0.84 eV the quasi- 
particle correction amounts to more than 2 eV. It can be 
expected that a treatment beyond the one-shot approach. 
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for example, by taking into account ofFdiagonal elements 
of the self-energy, by updating the Green function, or 
by including vertex corrections, will further increase the 
value and, thus, bring it even closer to the experimental 
value. As was shown in Ref. 1, already using LDA+J7 
instead of LDA as the mean-field starting point, which 
corrects the 2p-3d hybridization, gives an upward correc- 
tion of 0.2 eV in the resulting GW band gap. 

In conclusion, the band convergence is a serious issue 
in GW calculations and must be thoroughly dealt with. 
ZnO is an extreme case in this respect. We have reexam- 
ined the calculation of the ZnO band gap by Shih et al.^ 
and could confirm their main result: the GW band gap 
of ZnO shows a very slow convergence with respect to 
the number of states used to construct the polarization 
function and the correlation self-energy. Furthermore, 



when high-lying bands are included in the calculation, 
the linearization error of all-electron approaches becomes 
another important issue. As we have shown, it can be 
eliminated systematically within the FLAPW method by 
augmenting the basis in the MT spheres with LOs. In the 
case of ZnO this procedure yields a correction of about 
0.5 eV and brings the calculated band gap (2.99 eV) much 
closer to experiment than in previous studies. We believe 
that our study helps to clarify the contradiction between 
the pseudopotential results of Ref. [H and the previous 
results based on all-electron approaches!^— 
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